Fractional Fokker-Planck subdiffusion in alternating fields 
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The fractional Fokker-Planck equation for subdiffusion in time-dependent force fields is derived 
from the underlying continuous time random walk. Its limitations are discussed and it is then 
applied to the study of subdiffusion under the influence of a time-periodic rectangular force. As 
a main result, we show that such a force does not affect the universal scaling relation between 
the anomalous current and diffusion when applied to the biased dynamics: in the long time limit 
subdiffusion current and anomalous diffusion are immune to the driving. This is in sharp contrast 
with the unbiased case when the subdiffusion coefficient can be strongly enhanced, i.e. a zero- 
frequency response to a periodic driving is present. 

PACS numbers: 05.40.-a, 05.40.Fb, 05.60.-k, 02.50.Ey 
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I. INTRODUCTION 

The theoretical study of anomalously slow relaxation 
processes in time-dependent force fields constitutes a 
challenge of current research interest which is not free 
of ambiguity. It is known that there is no unique phys- 
ical mechanism responsible for the occurrence of subd- 
iffusion in condensed media One possible mecha- 
nism, which will be addressed in this work, corresponds 
to disordered glassy-like media consisting of trapping do- 
mains where the traveling particle can dwell for a random 
time with divergent mean value 0, d, 0|- The succes- 
sive residence times in traps are assumed to be mutually 
uncorrelated. Diffusion is nevertheless a non-Markovian 
(semi-Markovian) process exhibiting long (quasi-infinite) 
time correlations in the particle positions with a weak er- 
godicity breaking |[5j]. Mathematically this physical pic- 
ture can be described by a continuous time random walk 
(CTRW) model di] which in the continuous space limit 
leads to the fractional Fokker-Planck equation (FFPE) 
[1, 0j 9- This latter formulation is incomplete, as no 
non-Markovian master equation can define the underly- 
ing non-Markovian stochastic process [9]. However, the 
FFPE is very useful and has a tightly associated, com- 
plete description with (ordinary) Langevin equation in 
subordinated, random operational time p^. [Til [l^. 

The generalization of the CTRW and FFPE to time- 
dependent forces is a highly non-trivial matter since the 
force changes in the real and not in the operational time 
0, Also, how a field varying in time affects the dis- 
tribution of the residence times in the traps is not clear 
without specifying a concrete mechanism or some plausi- 
ble model, especially when the mean residence time does 
not exist [l^. The FFPE describing the dynamics in 
time-dependent force fields Fix, t) becomes ambiguous 
with a frequently (ab)used, ad hoc version [13, [13 which 
lacks clear theoretical grounds fl3|. The correct version 
of the FFPE for time-dependent fields was first given 
in Kefs. 0, HI]: differently from the FFPE for a time- 
independent force, in the case of a time-dependent field 



the fractional derivative does not stand in front of the 
Fokker-Planck operator but after it. As we explain with 
this work in more detail, such a FFPE can be justified be- 
yond the linear response approximation within a CTRW 
approach only for a special class of dichotomously fluctu- 
ating fields. The derivation of the FFPE for subdiffusion 
in such time-dependent fields is presented in Sec. [TTl 



In Sec, mil we apply the derived FFPE to study the in- 
fluence of time-periodic rectangular flelds on subdiffusive 
motion. Analytical solutions of the FFPE are conflrmed 
by stochastic Monte Carlo simulations of the underlying 
CTRW. In particular, we show with this work that the 
universal scaling relation between the biased anomalous 
diffusion and sub-current 0,0, 13, Hi is not affected by the 
periodic driving. Neither current nor diffusion are influ- 
enced asymptotically by the time-periodic field. This is 
in spite of the fact that unbiased subdiffusion of the stud- 
ied kind can be strongly enhanced in the time-periodic 
field [ll[ll]. 



II. DERIVATION OF THE FFPE FOR 
TIME-DEPENDENT FIELDS FROM THE 
UNDERLYING CTRW 



Since the FFPE does not define the underlying stochas- 
tic non-Markovian process, its generalization to include 
the influence of a time-dependent fleld should start from 
the underlying CTRW ^Q, [Tjl] . Following the general pic- 
ture of the CTRW we introduce a one-dimensional lat- 
tice {xi = iAx} with a lattice period Ax and i = 
0, ±1, ±2, . . . Let us first assume that there is no time- 
dependent field. After a random trapping time t a parti- 
cle at site i hops with probability qf to one of the nearest 
neighbor sites z ± 1; qf + q~ = 1. The random time r is 
extracted from a site-dependent residence time distribu- 
tion (RTD) ipiir). The corresponding generalized master 



2 



equation for populations Pi (t) reads [1^, [20, [2l| 



gf{t) = qf{t)vf. Then i^j(s) = vfs'^'" and we get 



(1) 



(3) 



The Laplace-transform of the kernel Kf{t) is related to 
the Laplace-transform of the RTD via Kf(s) = 

qf s'ipi{s)/[l — V'i(s)]- In the presence of a time- 
dependent field the kernels become generally functions 
of both instants of time and not only of their differ- 
ence, i.e. Kf{t — t') Kf{t,t'). One can relate 
Kf^{t,t') with the corresponding time-inhomogeneous 



^(rji), which are conditioned 
However, one always needs 



'[9tit')+97it')]dt']e.nd 
This yields in Eq. ((H) 



RTDs ipf{t + T,t) = i)\ 
on the entrance time t. 
a concrete an d p hysically meaningful model to pro- 
ceed further [1^. A simple example is a Marko- 
vian process with time-dependent rates g^{t), where 

V'f(r|i)-.9f(t+r)exp{-/;+^ 

Kf{t,t') = 2gHt)S{t-t') 
the standard master equation for a time-inhomogeneous 
Markovian process. How the non-exponential RTDs will 
be modified for a time-inhomogeneous process is gener- 
ally not clear In the present case, one can assume 
that the trapping occurs due to the existence of direc- 
tion(s) orthogonal to the x-coordinate. According to 
the modeling in Refs. an external field directed 

along X would not affect the motion in the orthogonal 
direction(s). However, it is not correct to think that the 
RTD in the trap will not be influenced by the field act- 
ing in the direction of a;, as it will change the rates (let 
us assume this simplest, tractable model) for moving left 
or right when escaping from the trap. Therefore, the 
RTD will generally be affected, see e.g. Ref. \2^. Ob- 
viously, the only situation when the RTD in the trap 
will not be changed is when the sum of the rates to 
escape from the trap, either left or right, is constant. 
In that case, only the probabilities q^{t) acquire addi- 
tional time-dependence and not the RTDs ipiir). This 
corresponds to the special class of dichotomously fluctu- 
ating force fields F{xi,t) = F{xi)£_{t), where ^{t) = ±1. 
Beyond this class, at most the linear response approxi- 
mation can work (l6j . Therefore, we restrict our treat- 
ment to the above class of fluctuating potentials. In this 
case, we can write Kf(t,t') = qf{t)Ki{t — t'), where 
Ki(s) = s'4'i{s)/[\ — ipi{s)]. Furthermore, we use the 
Mittag-Leffler distribution for the residence times 



(2) 



Here Ea{z) = I],T=o -^"/rC"" + 1) denotes the Mittag- 
Leffler function, a £ (0, 1) is the index of subdiffusion, 
and Vi = [gt{t) + gi{t)Y/"' the time scaling parameter; 



where the symbol o^t stands for the integro- 
differential operator of the Riemann-Liouville fractional 
derivative acting on a generic function of time x(t) as 



dt' ■ 



x(t') 



I d 



T{a) dt Jo {t - t'Y 



(4) 



r(a) is the gamma- function. In a time-dependent poten- 
tial U{x,t), one can set 

gfit) = iK^/Ax^)eM-m±i/2it)-UM} 

« {K^/Ax'')exp[±PFix„t)Ax/2], (5) 

so that the Boltzmann relation gt-i[i) / gi i^) = 
exp{/3[[/i_i(t) — Ui{t)\} is satisfied exactly and the time- 
independence of gf{t) +g^{t) = = const is also main- 
tained for small Ax and a sufficiently smooth potential. 
We have used here the notation /7i(t) = U{iAx,t) and 
Ui±i/2{t) = U{iAx±Ax/2,t); (3 = ksT is the inverse of 
temperature and Kq free fractional diffusion coefficient 
with dimension cm^s^". By passing to the continuous 
space limit Aa; as in Ref. f^, one finally obtains, 



d_ 

di 



P{x,t) 



a F{x,t) 

dx rja " dx^ 



oD] V(x,t).(6) 



In the latter equation r]a = il3Ka)~^ is the fractional 
friction coefficient. 

In the following we use Eq. ^ to study analyti- 
cally the subdiffusion in time-periodic rectangular fields. 
Our study is complemented by stochastic simulations of 
the underlying CTRW using the algorithm detailed in 
Ref. i. 



III. DRIVEN SUBDIFFUSION 

We consider a dichotomous modulation of a biased sub- 
diffusion where the absolute value of the bias is fixed but 
its direction flips periodically in time, i.e. 



Fit) = Fom 



(7) 



with 



_ / +1 foi' riTQ <t < {n + r)To 

\-l for {n + r)To < t < {n + 1)to ' 

Here tq is the period of the time-dependent force and 
n — 0,1,2... The quantity r G (0,1) determines the 
value of the average force: 



F = {Fit)U = Foi2r - 1) . 



(9) 
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FIG. 1: (Color online) Average particle position {x{t)) for var- 
ious values of the parameter r (average force F) and anoma- 
lous exponent a. Symbols represent the results from the nu- 
merical simulations of the CTRW obtained by averaging over 
10^ trajectories (for r — 0.5 over lO'' trajectories). Continu- 
ous lines represent the analytical solution (|14p of the FFPE 
The time-period of the rectangular force ((?]) is tq — 1 
and the fractional exponent a = 0.5; in numerical simulations 
-P<) / (??Q xA^q") = 1 is used. The value r — 0.9 corresponds to 
F = 0.8^0 and a = O.GFo; r = 0.8 corresponds to F = O.6F0 
and a — O.8F0; r = 0.7 corresponds to F — OAFo and 
a ^ O.92F0; r = 0.6 corresponds to F = Q.2Fo, a ^ 0.98f(); 
and r = 0.5 to F = 0, a = Fq. 



For r = 0.5 the average bias is zero and we recover the 
model investigated in Ref. Notice that the force 

F{t) can be decomposed in the following way: F{t) = 
F + F{t). The asymmetric driving, 



r) for UTo <t < {n + r)To . s 

for {n + r)ro < < < (n + 1)to ^ ' 



has a zero mean value, {F{t))ro — 0, and the driving 
root-mean-squared (rms) amplitude is c = {F'^{t)) 
2Fo-\/r(l — r). For a fixed average bias F, this yields 



1/2 



2r - 1 



(11) 



and therefore one can vary the ratio a/F between for 
r — 1 and 00 for r = 1/2-1- e, e 0. This offers the 
way to study the influence of an asymmetric, zero-mean 
driving F{t) with period tq and rms amplitude a on the 
subdiffusion under constant bias F. 

Let us begin by flnding the recurrence relation for the 
moments Assuming in Eq. ^ the force of the 

form ^ with ([8]), multiplying both sides of Eq. ([6]) by 
cc", and integrating over the x-coordinate one obtains, 

+ n(n-l)«:„o£>r"(-^""'W), (12) 
with subvelocity Va = -fo/^/a {n > 1). For n — 1 the last 



FIG. 2: (Color online) The analytical solution p4p for av- 
erage particle position {x{t)) obtained from the FFPE ((6)1 is 
presented for various values of the parameter r (average bias 
F) and anomalous exponent a. The time-period of the force 
is To = 1, however, in the long time limit the same asymptotic 
value is obtained for any value of tq. The relation between r 
and F and a is the same as in Fig. [1] 



term on the right hand side of Eq. (|T2|) is absent, 
dixit)) 



dt 



Via) 



(13) 



Equations (|T2)) and (|T3l) will be used to calculate the 
average particle position and the mean square displace- 
ment. 



A. Average particle position 

Upon integrating Eq. in time with ^(t) given by 
the solution for the average particle position reads: 



<t)) 



+ , NTo<t<iN + r)ro 



Vat' 



iN + r)To<t<iN + l)To 



(14) 



with 



XN = ixiO)) 



Tia + l) r(a-fl) 



^ [2(n + r)"-n"-(n+l)"] , (15) 
iN + r)"; (16) 



n=0 



, 2?;q,Tq 
Xn = xn + 



Tia + l) 

N counts the number of time periods passed. 

When the average bias is zero, i.e., r = 0.5, in the 
long time limit the mean particle position approaches 
the constant value 



ixi^)) ^Va,T^b{a)/T{a + l), 



(17) 
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with b{a) = Ejr=o[2(" + 1/2)" - - (n + 1)«]. The 
function b{a) changes monotonously from 6(0) = 1 to 
6(1) — 0. It describes the initial field phase effect which 
the system remembers forever when a G (0, 1) (see also 
Ref. [13]). This is one of the main differences between 
the anomalous motion in the absence of a force and in 
the presence of a time-dependent field with zero average 
value. 

In Fig.[l]the analytical solution p4)) for the mean par- 
ticle position {x{t)), obtained from the FFPE ((6|), is com- 
pared with the numerical solution of the CTRW for dif- 
ferent values of r, i.e. for different values of the average 
bias F. In Fig. [2] the solution (|14p is presented in the 
long time limit for various values of r and a. Figures [T] 
and[2]demonstrate that in the presence of an average bias 
the mean particle position grows as i". For all values of 
r, the asymptotic value of subvelocity corresponds to the 
averaged bias F = Fo(2r — 1), indicating that the periodic 
unbiased field does not affect the subdiffusion current for 
different values F and the field rms-amplitude a. 

Furthermore, the results depicted in Figs. [1] and 
clearly show the phenomenon of the "death of linear re- 
sponse" of the fractional kinetics to time-dependent fields 
in the limit t —f oo: the amplitude of the oscillations de- 
cays to zero as l/t^"", Eq. ^ (see also Refs. [l6l[23j). 
The amplitude of the oscillations is larger for larger val- 
ues of a and of tq. However, r(l -I- a){x{t)) /t°' reaches 
asymptotically the same value for any a and tq. 



B. Mean square displacement 

Let us now study the mean square displacement, de- 
fined as 



For n — 2 one obtains from Eq. (|12p . 



(18) 



dt 



2^n 



r(a) 




FIG. 3: (Color online) Scaled effective fractional diffusion co- 
efficient K,a^^ versus fractional exponent a for different driv- 
ing periods tq. The average bias is zero (r = 0.5). The an- 
alytical prediction (I23|l (continuous lines) is compared with 
the results (symbols) obtained from the numerical simula- 
tion of the CTRW by averaging over 10^ trajectories. For 
To > To ~ 8.818 [l3l the effective fractional diffusion coeffi- 
cient fi:i°^' {a) exhibits a maximum. 



time limit (see Appendix f^ . 



{x^{t)) 



2«2(2r 



1)2 



r(2Q!-i- 1 

2vl{2r-l)Si 
io^T{a + 1) 

C(2 + a 



r(a-t- 1) 
8^2 cos(a7r/2) 
7r2a;^r(a -I- 1) 



cos(nr27r) 



n=l 



(22) 



here C,{x) is the Riemann's zeta-function and 5*1 is a func- 
tion of a and r as given by Eq. (jA7| in Appendix [Al 

For r = 0.5 (average zero bias) the first and third term 
in Eq. are equal to zero. Furthermore, in the long 
time limit the average particle position (a;(oo)) is a finite 
constant. The asymptotic behavior of the mean square 
displacement is thus proportional to as in the force 
free case, however, characterized by an effective fractional 

Jeff) 



In order to find the analytical solution for the mean diffusion coefficient Ka instead of the free fractional 

square displacement, we use the Laplace-transform diffusion coefficient Kq, i.e. {5x'^{t)) = 2Ka t°- lT{l^a) 

method and the Fourier series expansion for S,{t) = for t ^ oo. The effective diffusion coefficient is |l3]. 
^{t + To) given wuth Eq. (gl). 



£.{t) ^ ^ /„exp(mwot). 



(20) 



4°^^ = -a + C(a + 2) 



with 



1 - 



1 



2"+2 



cos(a7r/2) 



(23) 



1 r" 

fn = — £_{t) exp{-mu}ot)dt 

To Jo 

= [1 — cxp(— mr27r)]/(m7r) (21) 

and loq = 27r/ro. Applying them to Eq. (I19p and assum- 
ing {x{0)) = and (a;2(0)) = we obtain that in the long 



The driving-induced part of the effective subdiffusion co- 
efficient is directly proportional to the square of driv- 
ing amplitude Fq and inversely proportional to lOq. For 
slowly oscillating force fields this leads to a profound ac- 
celeration of subdiffusion compared with the force free 
case: an optimal value of the fractional exponent a ex- 
ists, at which the driving-induced part of the effective 



5 



to 




FIG. 4: (Color online) The asymptotic scaling relation (|25|l . 
Symbols correspond to the numerical results obtained from 
the CTRW for different values of r (r / 0.5) and tq. Solid 
curve corresponds to the analytical result (|25p . Different val- 
ues of the parameter r correspond to different values of the 
bias F and the periodic field rms a, as described in Fig. [1] 
The field parameters do not influence the results within the 
statistical errors. 



fractional diffusion coefficient possesses a maximum (see 
Fig. ED. 

When r ^ 0.5 (finite average force) we obtain in the 
long time limit for {x{t))'^ (see Appendix [B|) . 



^a(2r-l)- 
r2(a + 1) 



^2" + 



2vl{2r-l)S, 
w^r(a + 1) 



(24) 



\Si is given by Eq. (|A7[) in Appendix [X] . Clearly, the 
leading term in Eq. (|24p corresponds to subvelocity in 
constant field F (averaged bias), i.e. the influence of 
periodic, unbiased driving F(t) dies asymptotically out, 
as illustrated in Fig. [21 

The results (HH) and ([M)) indicate that in the presence 
of a rectangular time-periodic force with a finite average 
value the general behavior of the mean square displace- 
ment is similar to the case of a constant force, i.e. the 
mean square displacement {Sx'^(t)) consists of terms pro- 
portional to and t^". In fact, for the leading term 
proportional to i^" in the mean square displacement one 
obtains the coefficient 



1 



r(2a + l) r2(a + l) 



This coefficient is the same as in the case of the subdifFu- 
sive motion under the influence of a constant force if the 
value of the constant force would be F. 

Furthermore, similarly to the case of a constant bias, 
the asymptotic scaling relation holds (r ^ 0.5, F 0) be- 
tween the mean square displacement and average particle 
position. In the limit t — > cx) the mean square displace- 
ment grows as {Sx^it)) (x t^" and 



A"o^ (x(t))2 ~ r(2a+ 1) 



lim 



1 



Fj^O. 



(25) 



It is illustrated by Fig. [4l where the analytical curve 
[Eq. (|25p ] is compared with the numerical results. The 
universality of relation ((25|) under the unbiased driving 
F{t) means that the biased diffusion is not affected by 
the driving. This is in a sharp contrast with the unbi- 
ased diffusion in Fig. [3l 



IV. CONCLUSION 

With this work we presented the derivation of the 
FFPE ([6]) for a special class of space- and time-dependent 
force fields from the underlying CTRW picture. Our 
derivation shows along with the corresponding discus- 
sion that it is difficult to justify this equation for time- 
dependent forces different from F{x,t) = F{x)^{t) with 
^{t) = ±1 beyond the linear response approximation. 
Using the FFPE ^ we demonstrated that the uni- 
versal scaling relation ((25)) for a biased subdiffusion is 
not affected by the additional action of a time-periodic 
zero-mean rectangular driving; neither is the asymptotic 
anomalous current nor the anomalous biased diffusion. 
We argue that this result is general and it is valid for 
other driving forms as well. This driving-immunity is due 
to the fact that the CTRW subdiffusion occurs in a ran- 
dom operational time which lacks mean value, whereas 
any physical field changes in the real, physical time. The 
CTRW-based subdiffusion fails to respond asymptoti- 
cally to such time-dependent fields while on its intrinsic 
random operational time scale any real, alternating field 
is acting infinitely fast and it makes effectively no influ- 
ence in a long run [l^, [3] , unless the rate of its change is 
precisely zero. This is the main reason for the observed 
anomalies. The remarkable enhancement of the unbi- 
ased subdiffusion within the CTRW framework by time- 
periodic rectangular fields is rather an exception than the 
rule. 
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APPENDIX A 

Using the property C{df{t)/dt} = sC{f{t)}-f{0) and 
assuming the initial conditions {x{t)) = and {x'^{t)) = 
we obtain from Eq. p9)) . 



C{{x^t))} 



m)}- (Ai) 
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Considering that 



here ({x) is the Riemann's zeta-function. Analogously, 



C{oDt "(xW)} = s'-"C{{xm = s-"C{d{x{t))/dt} , 
one obtains, 

fn 



71 — — OO 



^ ^s—iriLOQ 



Using Eq. (HH) with it follows, 



d{x{t)) \ 
dt J 



fn 



m— — OO 



[s — iLUo{n + to)]" 



Inserting (K2^ and (K3^ into Eq. (|XT|) we obtain. 



(A3) 



9 9 2 



s ^ — ' (s — mcjn)" 

n— — OO ^ ^ 



E 



/rr 



(A4) 



m— — OO 



Let us separate in Eq. (|A4[) the terms m = and n = 0, 



s 



ct+1 



/n 



n 2 °° 

(s — imcjn)" s ^ (s — inujn)" 

m— — OO ^ n— — OO ^ ^ 



E 



[s — iujoln + 771)]"= 



(A5) 



In the long time limit, i.e., in the limit s 0, in the 
double sum only terms with to = — n contribute, giving 
thus. 



?r) — — r 



m/0 



^ n — — OO ^ ' 

Let us compute the sums. Considering that 



(A6) 



E 

m— — OO 



fry 



{—imY 



E 



f-r 



and replacing here fm from (|2ip . one obtains. 



E 

m — — OO 



/n 



(— ito)" TT 



C(l + a) sin (Q;7r/2) 



E 

m— 1 



sin [(a — 4TOr)7r/2] 



TO 



l+a 



5i ; (A7) 



E 

n— — OO 



I /"I 4 , , , 

— — ^ cos(Q;7r/2) 



(—in)" TT^ 



C(2 + a)-^ 



cos(nr27r) 



(A8) 



Replacing these sums into Eq. (jA6p and considering that 
/o = 2r - 1 [see ([21]) together with dH])], we get. 



oQ + l 



+ , cos(a^/2) 



C(2 + a)-^ 



cos(7ir27r) 



2+Q 



(A9) 



Taking here the inverse Laplace transform one obtains 
the expression for {x'^{t)) in the long time limit [Eq. 



APPENDIX B 



Using Eq. ([13]), the quantity {x{t))'^ can be written in 
the following way. 



(xit))' 



2?;2 /■* /■*' 



rna) Jg Jo 



(Bl) 

Exploiting the property C{J^ f{t')dt'} = s-'^C{fit')} 
and denoting t' = t and = i', we can write, 



r2(a) 



2u; 



E " E 



r2(a)s ^ 



n— — OO m— — OO 



X <; / exp(imc^oi')rfi' J' ■ (B2) 



For a = 1 the latter equation gives. 



= V — 

^ ^' s ^ s- inojQ 

n— — OO 



E 



m— — OO 



In 



s — iLUo{n + to) 



(B3) 



Comparing this result with Eq. (|A4|) wc sec that 
£{(5x2(t))} = 2Ka/s'^, and thus {dx'^{t)) — 2^^^, as it 
should be for normal Brownian motion. 

For < a < 1 it is more convenient to proceed as 
follows. Let us calculate the Laplace trasform of {x{t)). 



7 



Considering that C{{x(t))} = s'^ C{d{x{t)) / dt} one ob- 
tains [see (|A3p ]. 



fn 



n— — oo 



s{s — inujo)^ 



n— — oo 
n/0 



fn 



s{s — intoo) 



In the hmit t —>■ oo, i.e. s — > 0, the latter equation 
becomes, 



Taking into account (|A7p we can write. 



(B6) 



Taking here the inverse Laplace transform we have for 
t oo, 



r(a + 1) 



OJr, 



C{{x{t))} 



E 



fn 



S ^ — ' ( — mCJn)" 

TL— — 00 



(B5) From here one obtains the result for {x{t))^ 
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